DNA methylation analysis

cat(paste(Sys.Date()))

## 2018-05-31

Load libraries

library(knitr)
library(data.table)
library(tidyverse)
library(minfi)
library(minfiData)
library(sva)
library(doParallel)
library(limma)
library(DMRcate)
library(IlluminaHumanMethylation450kanno.ilmn12.hg19) 
library(IlluminaHumanMethylation450kmanifest)
library(RColorBrewer)
library(ChIPpeakAnno)
library(EnsDb.Hsapiens.v75)
library(biomaRt)
library(pheatmap)
library(ggplot2)
library(gridExtra)

Helper functions

Normalization

## Run ComBat on a GenomicRatioSet to remove batch effects.
## In:
## - A GenomicRatioSet
## - An array with the same length as the number of experiments 
##   in the GenomicRatioSet, with the batch for each experiment.
## Out:
## - A new GenomicRatioSet, after correction of batch effects.
combat.on.grset <- function(GRset, batch){
    pheno <- pData(GRset)
    combat.model = model.matrix(~1, data=pheno)
    Ms_combat=ComBat(dat=getM(GRset), batch=batch, mod=combat.model)    
    
    tmp = RatioSet(M=Ms_combat, metadata=pheno)
    annotation(tmp) = annotation(MsetEx)
    GRset = mapToGenome(tmp)
    
    return(GRset)
}

# Use minfi quantile normalization, then filter probes and correct batch effects.
## Input:
## - A sample table
## - Remove probes with detection p-values above this cutoff (defualt is 0.01)
## - A table with non-specific probes
## - If probes with SNPS should be dropped (defualt is TRUE)
## - If we should run ComBat to correct for batch effects (defualt is TRUE)
## Output:
## - A GenomicRatioSet object, after quantile normalization and filtering of probes.
minfi.filter.and.normalize.by.tissue <- function(sample.table, detection.pval=0.01, drop.snps=TRUE, non.specific.probe.tab=NULL, run.combat=TRUE){
    RGSet <- read.metharray.exp(targets = sample.table)
    
    GRset.quantile <- preprocessQuantile(RGSet, fixOutliers = TRUE, removeBadSamples = TRUE, badSampleCutoff = 10.5, quantileNormalize = TRUE, stratified = TRUE, mergeManifest = FALSE, sex = NULL)

    # Filter on detection p-values  
    if(detection.pval<1){
        detP <- detectionP(RGSet)
        # ensure probes are in the same order in the mSetSq and detP objects
        detP <- detP[match(featureNames(GRset.quantile),rownames(detP)),]
        # remove any probes that have failed in one or more samples
        keep <- rowSums(detP < detection.pval) == ncol(GRset.quantile)
        GRset.quantile <- GRset.quantile[keep,] 
    }

    # Remove probes overlapping SNPs
    if(drop.snps){
        GRset.quantile  <- dropLociWithSnps(GRset.quantile , snps=c("SBE","CpG"), maf=0)
    }   
    
    # Remove non-specific probes
    if(!is.null(non.specific.probe.tab)){
        keep <- !(featureNames(GRset.quantile) %in% non.specific.probe.tab$TargetID)
        GRset.quantile <- GRset.quantile[keep,] 
    }
    
    if(run.combat){
        GRset.quantile <- combat.on.grset(GRset.quantile, batch=pData(GRset.quantile)$Sentrix_ID)
    }
    
    return(GRset.quantile)
}

QC plot functions

## Input:
## - A GenomicRatioSet object
## - Which tissue we are looking at (used as header)
## Output:
## - Density plot and density bean plot
minfi.qc.plot.by.tissue <- function(GRset,  sample.sheet, tissue="", sleep.column="Sleep"){
    densityPlot(getBeta(GRset), sampGroups = factor( sample.sheet[,sleep.column]), main=tissue)
    
    op <- par(mar=c(5,7,4,2)+0.1)
    densityBeanPlot(getBeta(GRset), sampGroups = factor(sample.sheet[,sleep.column]),
                                    sampNames=sample.sheet$Sample_Name, main=tissue)
    par(op)
}

Methods for DMRs

## Uses the DMRcate algorithm to find differentially methyltaed regions.
## In:
## - A GenomicRatioSet object
## - A string describing which model to use
## - A string describing the contrast to consider
## - FDR cutoff on differentially methylated probes (default is 0.05)
## - Fold change cutoff on differentially methylated regions (default is 0.02)
## - P-value (Stouffer) cutoff on differentially methylated regions (default is 0.05)
## - lamda: bandwidth to gaussian kernel smoothing function (default is 1000)
## - C: Scaling factor for bandwidth
## Out:
## - A GRanges object with info on all DMRs.
dmrcate.dmrs <- function(GRset, sample.sheet, use.model, contrast, cpg.fdr=0.05, dmr.fc=0.02, dmr.fdr=0.05, lambda=1000, C=2){
    designMatrix <- model.matrix(as.formula(use.model), data=sample.sheet) 

        
    contMatrix <- makeContrasts(contrasts=contrast, levels=designMatrix)
    
    myAnnotation <- cpg.annotate(GRset, datatype = "array", analysis.type="differential", design=designMatrix, arraytype="450K", contrasts = TRUE, cont.matrix = contMatrix, coef=contrast, fdr = cpg.fdr)
    
    
    ## If no DMRs are found, dmrcate throws an error. We need to catch this so the programs doesn't crash.
    ## Instead we now return an empty GRanges object.
    DMRs <- tryCatch(
        dmrcate(myAnnotation, lambda=lambda, C=C), 
        error = function(e)
            { ## Of no DMRs found, return NULL
            print(e$message)
            return(GRanges())
        }
    )
    
    if(!(class(DMRs) == "GRanges")){
        results.ranges <- extractRanges(DMRs, genome = "hg19")
        results.ranges <- results.ranges[abs(results.ranges$meanbetafc) > dmr.fc & results.ranges$Stouffer < dmr.fdr,]
        
        return(results.ranges)
    }
    else{
        return(DMRs)
    }
}

## Uses limma algorithm to find differentially methyltaed CpGc/probes.
## In:
## - A GenomicRatioSet object
## - A string describing which model to use
## - A string describing the contrast to consider
## - The number of most significant CpGs to return
## Out:
## - A a table with the most significant CpGs, and the following columns:
##  ID, CHR, pos, logFC, AveExpr, t, P.Value, adj.P.Val, B
diff.cpgs <- function(GRset, sample.sheet, use.model, contrast, nr.cpgs=20){
    designMatrix <- model.matrix(as.formula(use.model), data=sample.sheet) 
    contMatrix <- makeContrasts(contrasts=contrast, levels=designMatrix)
    
    ## Annotate CpGs.
    cpg.fdr=1
    myAnnotation <- cpg.annotate(GRset, datatype = "array", analysis.type="differential", design=designMatrix, arraytype="450K", contrasts = TRUE, cont.matrix = contMatrix, coef=contrast, fdr = cpg.fdr)
    
    myAnnotation.tab <- as.data.frame(lapply(myAnnotation,head, n=1e6))
    rownames(myAnnotation.tab) <- myAnnotation.tab$ID

    ## Run limma on each CpG. 
    fit <- lmFit(getBeta(GRset), designMatrix)
    fit2  <- contrasts.fit(fit, contMatrix)
    fit2  <- eBayes(fit2)
    
    ## Format output
    tt <- topTable(fit2,coef=1, number=nr.cpgs)
    tt <- cbind(myAnnotation.tab[rownames(tt),c("CHR", "pos")], tt)
    tt$logFC <- signif(tt$logFC, 2)
    tt$AveExpr <- signif(tt$AveExpr, 2)
    tt$P.Value <- signif(tt$P.Value, 2)
    tt$adj.P.Val <- signif(tt$adj.P.Val, 2)
    tt$B <- signif(tt$B, 2)
    return(tt)
}

Functions to plot individual probes and DMRs

## Get probes included in a DMRs found by bumphunter
## In:
## - Which program was used to get DMRs ("bumphunter" or "dmrcate")
## - DMRs (bumphunter output, bump object or DMRcate output, GRranges object)
## - The GenomicRatioSet on whic bumphunter was run
## - An index decribing wgich DMR we are interested in
## - Direction of change (1 is more methylation, -1 is less methyltaion)
## Out:
## - An array with probe ids
get.dmr.probes <- function(method, dmrs, GRset, which.dmr, change.dir){
    dmr <- NULL
    if(method == "bumphunter"){
        dmr.tab <- dmrs$table[dmrs$table$value * change.dir > 0,]
        
        chr <- dmr.tab[which.dmr,"chr"]
        start <- dmr.tab[which.dmr,"start"]
        end <- dmr.tab[which.dmr,"end"]
        dmr <- GRanges(seqnames = chr, ranges = IRanges(start = start, end = end))
    }
    if(method == "dmrcate"){
        dmr <- dmrs[which.dmr,]
    }   
    
    probeinfo <- makeGRangesFromDataFrame(minfi::getAnnotation(GRset), start.field = "pos", end.field = "pos", keep.extra.columns=TRUE)
    GRanges.dmr <- probeinfo[from(GenomicRanges::findOverlaps(probeinfo, dmr, ignore.strand=TRUE)),]
    return( names(GRanges.dmr) )
}


## Get mean beta values for all probes in a differentially methylated region (DMR).
## Input:
## - The differentially methyltaed regions, as a GRanges object
## - The normalized methylation data, as a GenomicRatioSet object
## - Information on all probes, as a GRanges object created by makeGRangesFromDataFrame
## - Which DMR to get probes for, an intege index.
## Out:
## - An array mean beta values for a DMR, over all samples.
get.dmr.betas <- function(dmrs, GRset, probeinfo, which.dmr ){
    get.dmr.probes <- function(dmrs, GRset, which.dmr, probeinfo){
        dmr <- dmrs[which.dmr,]
        
        GRanges.dmr <- probeinfo[from(GenomicRanges::findOverlaps(probeinfo, dmr, ignore.strand=TRUE)),]
        return( names(GRanges.dmr) )
    }
    
    betas <- getBeta(GRset)
    dmr.probes <- get.dmr.probes(dmrs, GRset, which.dmr, probeinfo)
    probe.betas <- betas[dmr.probes,]
    dmr.betas <- apply(probe.betas, 2, mean)
    
    return(dmr.betas)
}


## Get differences in beta values for all probes in a differentially methylated region (DMR).
## Input:
## - The differentially methyltaed regions, as a GRanges object
## - The normalized methylation data, as a GenomicRatioSet object. ASSUMES THAT THE SAMPLES ARE
##   ORDERED AS FOLLOWS: subject1 control, subject1 treat, subject2 control, subject2 treat, ...
## - Information on all probes, as a GRanges object created by makeGRangesFromDataFrame
## - Which DMR to get probes for, an intege index.
## Out:
## - An with mean beta values differences for a DMR, for all subjects in the study.
get.dmr.betas.change <- function(dmrs, GRset, probeinfo, which.dmr ){
    get.dmr.probes <- function(dmrs, GRset, which.dmr, probeinfo){
        dmr <- dmrs[which.dmr,]
        
        GRanges.dmr <- probeinfo[from(GenomicRanges::findOverlaps(probeinfo, dmr, ignore.strand=TRUE)),]
        return( names(GRanges.dmr) )
    }
    
    betas <- getBeta(GRset)
    dmr.probes <- get.dmr.probes(dmrs, GRset, which.dmr, probeinfo)
    probe.betas <- betas[dmr.probes,]
    probe.beta.per.subject <- do.call("cbind", lapply(seq(from=1, to=30, by=2), function(i){ probe.betas[,i+1] - probe.betas[,i] }))
    dmr.betas.per.subject <- apply(probe.beta.per.subject, 2, mean)
    
    return(dmr.betas.per.subject)
}


## Plots differnces in methylation levels between conditions, for a singe probe of a DMR.
## Input:
## - A GenomicRatioSet
## - An id to plot, either a name of a probe in the GRset, or a name of a DMR in the GRanges object.
## - A Granges object describing DRMs. (Optional, deafualt is NULL in which case the id is expected to be a probe id.)
## - A title to put in the plot. (Default is NULL, in which case the id is used as a title).
## - Which samples correspond to the sleep condition (default 1,3,5,..,29)
## - Which samples correspond to the wake condition (default 2,4,6,..,30)
plot.methylation.difference <- function(GRset, select.id, dmrs=NULL, main=NULL, sleep=seq(from=1, to=30, by=2), wake=seq(from=2, to=30, by=2)){
    betas <- getBeta(GRset)
    selected.betas <- NULL
    
    # If dmrs data is provided, assume we want to plot average over a dmr.
    if(!is.null(dmrs)){
        dmr.probes <- get.dmr.probes(method="dmrcate" , dmrs=dmrs, GRset=GRset, which.dmr=select.id)
        probe.data <- betas[dmr.probes,]
        selected.betas <- apply(probe.data,2,mean)
    }
    else{
        selected.betas <- betas[select.id,]
    }
    
    if(is.null(main)){
        main <- select.id
    }
    
    plot.data <- data.frame(sleep=pheno.data$Sleep[c(sleep,wake)], subject=pheno.data$Subject[c(sleep,wake)], beta=c(selected.betas[sleep], selected.betas[wake]))
    
    p <- ggplot(plot.data ) +
        theme_classic() +
        geom_violin(aes(x = sleep, y = beta, fill=sleep), trim=FALSE) +
        geom_line(aes(x = sleep, y = beta, group = subject), size=0.2) +
        ggtitle(main) +
        xlab("")

    return(p)
}

Functions for pathway analysis etc.

## Given a set of genomic regions, resturns genes with TSS close to those regions.
## Input:
## - A set of genomic regions, GRanges
## - Max distance between TSS and regions (defualt=5000bp)
## - Annotation. GRanges object with 
## Output:
## - An array with (Ensembl) gene ids, of genes close to the given regions.
## - An array with rownumbers of genomic regions that are close to the given Ensembl genes.
annotate.regions.to.tss <- function(regions, maxgap=5000, annoData){
    overlaps.anno <- annotatePeakInBatch(regions, AnnotationData=annoData, featureType="TSS", FeatureLocForDistance="TSS", select="all", output="nearestLocation", maxgap=maxgap)
    
    overlaps.anno$ensg <- annoData$gene_id[match(overlaps.anno$feature, annoData$tx_id)]
    close.to.tss.indx <- abs(overlaps.anno$distancetoFeature)<= maxgap
    
    close.to.regions.ensg <- unique(overlaps.anno$ensg[close.to.tss.indx])
    close.to.genes.region <- unique(overlaps.anno$peak[close.to.tss.indx])

    peak.gene.tab <- data.frame(region=overlaps.anno$peak, ensg=overlaps.anno$ensg)
    peak.gene.tab <- peak.gene.tab[close.to.tss.indx,] ## Only keep rows where probe is close to tss
    peak.gene.tab <- peak.gene.tab[order(peak.gene.tab$region,peak.gene.tab$ensg),]
    
    ## Create list mapping ensg ids <- probe ids
    peak.gene.map <- split(x=as.character(peak.gene.tab$region), f=as.character(peak.gene.tab$ensg))

    return(list(close.to.regions.ensg, close.to.genes.region, peak.gene.map))
}


## Convert an array of Ensembl gene ids to Entrez ids. Genes with missing ids are removed.
## Input:
## - An array with ensembl ids
## - A table with columns ensembl_gene_id and entrezgene.
## Output: An array with Entrez gene ids
ensembl.set.to.entrez <- function(ensg.ids, gene.annot.tab){
    entrez.ids <- gene.annot.tab$entrezgene[match(ensg.ids, gene.annot.tab$ensembl_gene_id)]
    entrez.ids <- entrez.ids[!is.na(entrez.ids)]
    return(entrez.ids)
}


## Get enriched GO and KEGG terms, given a set of differentially methylated regions (DMRs).
## Input:
## - A set of DMRs (GRanges)
## - An array with (Ensembl) ids of all genes that have probes on the 450K array in the promoters.
## - A table with columns ensembl_gene_id and entrezgene.
## - Max Fisher p-value for enriched terms.
## - Minimum nr of genes for enriched terms.
## Output:
## - A list of enriched terms for: GO/BP, GO/MF, GO/CC and KEGG.
dmr.go.kegg <- function(dmrs, close.to.probe.ensembl, gene.annot.tab, annoData, fisher.cutoff=0.01, min.genes = 3){
    get.top.terms <- function(enrich.obj, fisher.cutoff, min.genes){
        enrich.obj[enrich.obj$P.DE< fisher.cutoff & enrich.obj$DE >= min.genes,]
    }
    
    dmr.genes <- annotate.regions.to.tss(dmrs, maxgap=5000, annoData)[[1]]
    
    close.to.probe.entrez <- ensembl.set.to.entrez(close.to.probe.ensembl, gene.annot.tab)
    dmr.genes.entrez <- ensembl.set.to.entrez(dmr.genes, gene.annot.tab) 
    
    go.fisher <- goana(dmr.genes.entrez, universe = close.to.probe.entrez, species = "Hs")
    kegg.fisher <- kegga(dmr.genes.entrez, universe = close.to.probe.entrez, species = "Hs", FDR = 0.05)
    
    out.list <- list(
        BP=get.top.terms(limma::topGO(go.fisher, number=100000, truncate.term=40, ontology = "BP"), fisher.cutoff, min.genes),
        MF=get.top.terms(limma::topGO(go.fisher, number=100000, truncate.term=40, ontology = "MF"), fisher.cutoff, min.genes),
        CC=get.top.terms(limma::topGO(go.fisher, number=100000, truncate.term=40, ontology = "CC"), fisher.cutoff, min.genes),
        KEGG=get.top.terms(limma::topKEGG(kegg.fisher, number=100000), fisher.cutoff, min.genes)
    )
    
    return(out.list)
}

Other functions

## Modified version of the function ChIPpeakAnno::assignChromosomeRegion.
## Instread of returning percentage and jacard index of regions annotated to different
## genomic features, an array with the feature for each input region is returned.
## Input:
## - the genomic regions to annoate (GRanges object)
## - upstream cutoff for promoters (default 1000)
## - downstream cutoff for promoters (default 1000)
## - Precedence of genomic features, in case a region overlaps several different fetures.
##   Has to be any combination of "Exon", "Intron", "fiveUTR", "threeUTR", "Promoter", 
##   and "immediateDownstream".
## - A txdb with the genomic features.
## Output: 
## - An array with the type of genomic feature for each input region.
## See ?assignChromosomeRegion for more info.
assignChromosomeRegion2 <- function (peaks.RD, proximal.promoter.cutoff = 1000L, 
                    immediate.downstream.cutoff = 1000L, precedence, TxDb) 
{
    if (!inherits(TxDb, "TxDb")) 
        stop("TxDb must be an object of TxDb, \n                     try\n?TxDb\tto see more info.")
    if (!inherits(peaks.RD, c("RangedData", "GRanges"))) 
        stop("peaks.RD must be a GRanges object.")
    if (!is.null(precedence)) {
            if (!all(precedence %in% c("Exon", "Intron", "fiveUTR", 
                                                                 "threeUTR", "Promoter", "immediateDownstream"))) 
            stop("precedence must be a combination of \nExons, Introns, fiveUTRs, threeUTRs, \nPromoters, immediateDownstream")
        }
    if (inherits(peaks.RD, "RangedData")){ 
        peaks.RD <- as(peaks.RD, "GRanges")
    }
    exons <- exons(TxDb, columns = NULL)
    introns <- unique(unlist(intronsByTranscript(TxDb)))
    fiveUTRs <- unique(unlist(fiveUTRsByTranscript(TxDb)))
    threeUTRs <- unique(unlist(threeUTRsByTranscript(TxDb)))
    transcripts <- unique(transcripts(TxDb, columns = NULL))
    options(warn = -1)
    try({
        promoters <- unique(promoters(TxDb, upstream = proximal.promoter.cutoff, 
                                                                    downstream = 0))
        immediateDownstream <- unique(flank(transcripts, 
                                                                                width = immediate.downstream.cutoff, start = FALSE, 
                                                                                use.names = FALSE))
    })
    microRNAs <- tryCatch(microRNAs(TxDb), error = function(e) return(NULL))
    tRNAs <- tryCatch(tRNAs(TxDb), error = function(e) return(NULL))
    options(warn = 0)
    annotation <- list(exons, introns, fiveUTRs, threeUTRs, promoters, immediateDownstream)
    if (!is.null(microRNAs)) 
        annotation <- c(annotation, microRNAs = microRNAs)
    if (!is.null(tRNAs)) 
        annotation <- c(annotation, tRNAs = tRNAs)
    annotation <- lapply(annotation, function(.anno) {
        mcols(.anno) <- NULL
        .anno
    })
    names(annotation)[1:6] <- c("Exon", "Intron", "fiveUTR", 
                                                            "threeUTR", "Promoter", "immediateDownstream")
    formatSeqnames <- function(gr) {
        seqlevels(gr)[grepl("^(\\d+|MT|M|X|Y)$", seqlevels(gr))] <- paste("chr", 
                                                                                                                                            seqlevels(gr)[grepl("^(\\d+|MT|M|X|Y)$", seqlevels(gr))], 
                                                                                                                                            sep = "")
        seqlevels(gr)[seqlevels(gr) == "chrMT"] <- "chrM"
        trim(gr)
    }
    peaks.RD <- formatSeqnames(peaks.RD)
    peaks.RD <- unique(peaks.RD)
    annotation <- lapply(annotation, formatSeqnames)
    annotation <- GRangesList(annotation)
    newAnno <- c(unlist(annotation))
    newAnno.rd <- GenomicRanges::reduce(trim(newAnno))
    Intergenic.Region <- gaps(newAnno.rd, end = seqlengths(TxDb))
    Intergenic.Region <- Intergenic.Region[strand(Intergenic.Region) != "*"]
    if (!all(seqlevels(peaks.RD) %in% seqlevels(newAnno))) {
        warning("peaks.RD has sequence levels not in TxDb.")
        sharedlevels <- intersect(seqlevels(newAnno), seqlevels(peaks.RD))
        peaks.RD <- keepSeqlevels(peaks.RD, sharedlevels)
    }
    mcols(peaks.RD) <- NULL
    annotation <- annotation[unique(c(precedence, names(annotation)))]
    names(Intergenic.Region) <- NULL
    annotation$Intergenic.Region <- Intergenic.Region
    anno.names <- names(annotation)
    
    ol.anno <- findOverlaps(peaks.RD, annotation)
    ol.anno.splited <- split(queryHits(ol.anno), anno.names[subjectHits(ol.anno)])

    ol.anno <- as.data.frame(ol.anno)
    ol.anno.splited <- split(ol.anno, ol.anno[, 2])
    hasAnnoHits <- do.call(rbind, ol.anno.splited[names(ol.anno.splited) != as.character(length(annotation))])
    hasAnnoHits <- unique(hasAnnoHits[, 1])
    ol.anno <- ol.anno[!(ol.anno[, 2] == length(annotation) & (ol.anno[, 1] %in% hasAnnoHits)), ]
    ol.anno <- ol.anno[!duplicated(ol.anno[, 1]), ]
    out.data <- names(annotation)[ol.anno$subjectHits] 
}

## Annotate DMRs, with genomic region and (possible) promoter.
## Input:
## - DMRs, as GRanges object
## - Max distance to transcription start site (default 5000 bp)
## - Annotation data as ??
## - Annotation data as txdb
## Output:
## - The same GRanges object as input, but with two additional columns: 
##   "gene_tss" (if the DMR is close to any TSS) and "genomic _feature"
##   ("near TSS", "exon", "intron", "Intergenic.Region" etc.)
annot.dmrs <- function(dmrs, maxgap=5000, annoData, annoData.TXDB){
    ## Find transcription start sites nearby the DMRs
    overlaps.anno <- annotatePeakInBatch(dmrs, AnnotationData=annoData, featureType="TSS", FeatureLocForDistance="TSS", select="all", output="nearestLocation", maxgap=maxgap)
    overlaps.anno$gene_name <- annoData$gene_name[match(overlaps.anno$feature, names(annoData))]
    overlaps.anno <- overlaps.anno[abs(overlaps.anno$distancetoFeature)<= maxgap, ] ## Only keep nearby TSS
    
    ## For each DMR, get names of genes with nearby TSS (if any)
    nearest.gene <- sapply(names(dmrs),
                                                 function(x){
                                                    peak.indx <- overlaps.anno$peak == x & base::grepl("ENST",overlaps.anno$feature)
                                                    paste(unique(overlaps.anno$gene_name[peak.indx]), collapse=", ")
                                                    })
    dmrs$gene_tss <- nearest.gene
    
    assigned.chromosome.regions <- assignChromosomeRegion2(dmrs, precedence=c("Promoter", "immediateDownstream","fiveUTR", "threeUTR", "Exon", "Intron"), proximal.promoter.cutoff=maxgap, immediate.downstream.cutoff=maxgap, TxDb=annoData.TXDB)
    
    dmrs$genomic_feature <- assigned.chromosome.regions
    
    ## Make sure the the tss annotations from annotatePeakInBatch and the assigned chromatin regions agree.
    gene.tss.overlap <- dmrs$gene_tss != ""
    dmrs$genomic_feature[gene.tss.overlap] <- "near TSS"

    return(dmrs)
}



## Summarize probe data by gene.
## Input:
## - The GenomicRatioSet with the probe data.
## - A list mapping gene id -> probe ids
## - Which function to use for summarizing. Default is mean.
## Output:
## - A matrix with gene as rows, samples as column, with summarized M-values of all probes annotated 
## to the gene in question. E.g. if the mean function is used, mean M-values over all probes for a 
## gene are returned.
mvals.per.gene <- function(GRset, probe.gene.map, use.function=mean){
    filtered.map <- lapply(probe.gene.map, 
                 function(x){
                    intersect(x, rownames(GRset))
                 })
    
    mvals <- getM(GRset)
    out.list <- lapply(names(filtered.map),
                                        function(ensg){
                                            probe.data <- mvals[filtered.map[[ensg]], ]
                                            
                                            if(length(nrow(probe.data)) ==0){
                                                return(NULL)
                                            }
                                            
                                            return( apply(probe.data, 2, use.function, na.rm=TRUE) )
                                        })
    names(out.list) <- names(filtered.map) 
    out.tab <- do.call(rbind, out.list)
    return(out.tab)
    }



## Export probe data as gzipped csv files with M-values.
## Input:
## - The GenomicRatioSet to export.
## - Which probes to export.
## - The name of the output file.
## Output:
## - Writes the M values from the given GenomicRatioSet to the given file.
write.m.vals <- function(GRset, use.probes, out.file){
    out.data <- getM(GRset)[rownames(GRset) %in% use.probes, ]
    out.data <- round(out.data,4)
    out.file.gz <- gzfile(out.file)
    write.csv(out.data , file=out.file.gz)
}

Load data and pre-process

# setwd()
set.seed(1977) ## set seed so that methods using randomization (bumphunter and SAM) always give the same results.

#########################################
## Minfi load and pre-process data

## Read info on non-specific probes
non.specific.probe.tab <- fread('https://raw.githubusercontent.com/sirselim/illumina450k_filtering/master/48639-non-specific-probes-Illumina450k.csv', header=TRUE)

## Read table with info about all arrays
targets <- read.csv(file="cedernaes2018_dnamethylation_samplesheet.csv")
rownames(targets) <- paste(targets$Sentrix_ID, targets$Sentrix_Position, sep="_")
targets.a <- targets[targets$Tissue=="A",]
targets.m <- targets[targets$Tissue=="M",]

## Normalize. Since there are large differences between adipose and muscle, use quantile normalization on each tissue separately.
minfi.norm.a <- minfi.filter.and.normalize.by.tissue (targets.a, non.specific.probe.tab=non.specific.probe.tab)

## Warning in .getSex(CN = CN, xIndex = xIndex, yIndex = yIndex, cutoff
## = cutoff): An inconsistency was encountered while determining sex. One
## possibility is that only one sex is present. We recommend further checks,
## for example with the plotSex function.

## Standardizing Data across genes

minfi.norm.m <- minfi.filter.and.normalize.by.tissue (targets.m, non.specific.probe.tab=non.specific.probe.tab)

## Warning in .getSex(CN = CN, xIndex = xIndex, yIndex = yIndex, cutoff
## = cutoff): An inconsistency was encountered while determining sex. One
## possibility is that only one sex is present. We recommend further checks,
## for example with the plotSex function.

## Standardizing Data across genes

QC plots

minfi.qc.plot.by.tissue(minfi.norm.a, sample.sheet=targets.a, tissue="Adipose (minfi)")

minfi.qc.plot.by.tissue(minfi.norm.m, sample.sheet=targets.m, tissue="Muscle (minfi)")

Find DMRs

## Get DMRs with dmrcate.
dmrcate.minfi.dmrs.a <- dmrcate.dmrs(minfi.norm.a, sample.sheet=targets.a, use.model="~ 0 + Sleep + Subject", contrast="Sleepwake-Sleepsleep", dmr.fc=0.02)
dmrcate.minfi.dmrs.m <- dmrcate.dmrs(minfi.norm.m, sample.sheet=targets.m, use.model="~ 0 + Sleep + Subject", contrast="Sleepwake-Sleepsleep", dmr.fc=0.02)

## [1] "The FDR you specified in cpg.annotate() returned no significant\n  CpGs, hence there are no DMRs.  Try specifying a value of\n  'pcutoff' in dmrcate() and/or increasing 'fdr' in\n  cpg.annotate()."

dmrcate.minfi.dmrs.nofc.a <- dmrcate.dmrs(minfi.norm.a, sample.sheet=targets.a, use.model="~ 0 + Sleep + Subject", contrast="Sleepwake-Sleepsleep", dmr.fc=0)
dmrcate.minfi.dmrs.nofc.m <- dmrcate.dmrs(minfi.norm.m, sample.sheet=targets.m, use.model="~ 0 + Sleep + Subject", contrast="Sleepwake-Sleepsleep", dmr.fc=0)

## [1] "The FDR you specified in cpg.annotate() returned no significant\n  CpGs, hence there are no DMRs.  Try specifying a value of\n  'pcutoff' in dmrcate() and/or increasing 'fdr' in\n  cpg.annotate()."

Find Cpgs

There were no signficant DMRs in muscle. Therefore we also look at the the most signficant individual probes/CpGs. None of these are very signficant either (from the FDR adjusted p-value P.adj):

diff.cpgs(GRset=minfi.norm.m, sample.sheet=targets.m, use.model="~ 0 + Sleep + Subject", contrast="Sleepwake-Sleepsleep", nr.cpgs=10)

##              CHR       pos  logFC AveExpr         t P.Value adj.P.Val   B
## cg02934600  chr1 107600376  0.028    0.34  6.774544 3.9e-06         1 4.0
## cg06070414  chr7  95226412  0.044    0.25  6.493057 6.6e-06         1 3.4
## cg00722081 chr10  61147668 -0.028    0.63 -6.327758 9.0e-06         1 3.1
## cg26493306  chr5 134463458  0.023    0.92  6.275672 9.9e-06         1 3.0
## cg04957146 chr17   2551652  0.021    0.87  6.164014 1.2e-05         1 2.8
## cg01051642  chr1   2517410  0.026    0.12  6.000548 1.7e-05         1 2.5
## cg09250087  chr2 161127564  0.029    0.69  5.862363 2.2e-05         1 2.2
## cg20292318 chr17  77141853 -0.025    0.82 -5.494034 4.5e-05         1 1.5
## cg10747185  chr3 176862267  0.038    0.84  5.493864 4.5e-05         1 1.5
## cg03063453 chr19  52693134  0.023    0.13  5.490893 4.5e-05         1 1.5

Analyze DMRs

Get numbers of DMRs

pval.cutoff <- 0.05
length(which(dmrcate.minfi.dmrs.a$meanbetafc > 0))

## [1] 62

length(which(dmrcate.minfi.dmrs.a$meanbetafc < 0))

## [1] 31

length(which(dmrcate.minfi.dmrs.m$meanbetafc > 0))

## [1] 0

length(which(dmrcate.minfi.dmrs.m$meanbetafc < 0))

## [1] 0

length(which(dmrcate.minfi.dmrs.nofc.a$meanbetafc > 0))

## [1] 92

length(which(dmrcate.minfi.dmrs.nofc.a$meanbetafc < 0))

## [1] 56

length(which(dmrcate.minfi.dmrs.nofc.m$meanbetafc > 0))

## [1] 0

length(which(dmrcate.minfi.dmrs.nofc.m$meanbetafc < 0))

## [1] 0

maxgap <- 5000
annoData <- toGRanges(EnsDb.Hsapiens.v75, feature="transcript")
annoData.TXDB <- makeTxDbFromGFF("/Users/orzechoj/projects/benedict2016/data/annot/Homo_sapiens.GRCh37.87.gff3.gz", format="gff3")

## Import genomic features from the file as a GRanges object ... OK
## Prepare the 'metadata' data frame ... OK
## Make the TxDb object ... OK

dmrcate.minfi.dmrs.a <- annot.dmrs(dmrcate.minfi.dmrs.a, maxgap=5000, annoData, annoData.TXDB)
dmrcate.minfi.dmrs.nofc.a <- annot.dmrs(dmrcate.minfi.dmrs.nofc.a, maxgap=5000, annoData, annoData.TXDB)

pheno.data <- targets.a
rownames(pheno.data) <- paste(pheno.data$Sentrix_ID, pheno.data$Sentrix_Position, sep="_")
probeinfo <- makeGRangesFromDataFrame(minfi::getAnnotation(minfi.norm.a), start.field = "pos", end.field = "pos", keep.extra.columns=TRUE)

## Plot all samples, DMRs with any fold change 
plot.data <- t(sapply(1:length(dmrcate.minfi.dmrs.nofc.a), function(i){get.dmr.betas(dmrcate.minfi.dmrs.nofc.a, minfi.norm.a, probeinfo, i)}))
plot.data <- plot.data[,c(seq(from=1, to=ncol(plot.data), by=2),seq(from=2, to=ncol(plot.data), by=2))]
pheatmap(plot.data, cluster_rows=TRUE, cluster_cols=FALSE, show_colnames=FALSE, scale="row", 
                 annotation_col=pheno.data[,c("Sleep","Subject")], main="Adipose DMRs (average Beta values)", fontsize_row=5)

Plot indivudial DMRs

## Which DMRs to print
plot.names <- c("CD36", "TSPAN", "GNAS", "INS", "GFI1", "AKR1CL1", "TNXB", "TRIM2", "HOXA2","FOXP2")
plot.dmrs <- names(dmrcate.minfi.dmrs.nofc.a[unlist(lapply(plot.names, function(x){grep(x,dmrcate.minfi.dmrs.nofc.a$gene_tss)})),] )
names(plot.dmrs) <- plot.names

## Make violin plots
p <- lapply(1:length(plot.dmrs),
                        function(i){
                            plot.methylation.difference(GRset=minfi.norm.a, select.id=plot.dmrs[i], dmrs=dmrcate.minfi.dmrs.nofc.a, main=names(plot.dmrs)[i])
                            })

ml <- marrangeGrob(p, nrow=5, ncol=2)
print(p)

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## Fetch gene info.
ensembl = useEnsembl(biomart="ensembl",GRCh=37, dataset="hsapiens_gene_ensembl")
gene.annot.tab <- getBM(attributes=c('ensembl_gene_id','gene_biotype','entrezgene', 'external_gene_name'), mart = ensembl)

## Clean up data: if there are several entries with same ensembl_gene_id, use the first
gene.annot.tab %>%
    group_by(ensembl_gene_id) %>%
    summarise_each(funs(paste(sort(.), collapse=","))) %>%
    extract(external_gene_name, "external_gene_name", "([^,]+),?.*") %>%
    extract(gene_biotype, "gene_biotype", "([^,]+),?.*") %>%
    extract(entrezgene, "entrezgene", "([^,]+),?.*") -> gene.annot.tab

## `summarise_each()` is deprecated.
## Use `summarise_all()`, `summarise_at()` or `summarise_if()` instead.
## To map `funs` over all variables, use `summarise_all()`

## Create genomic ranges object from info about 450K probes 
ann450k = minfi::getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19)
cpgData <- GRanges(seqnames=Rle(ann450k$chr),
                                     ranges=IRanges(start=ann450k$pos, end=ann450k$pos),
                                     strand=Rle(rep("*",nrow(ann450k))))
names(cpgData) <- rownames(ann450k)

## Get all genes that have 450K probes in promoters
probe.gene.overlaps <- annotate.regions.to.tss(cpgData, maxgap=5000, annoData)
close.to.probe.ensembl <- probe.gene.overlaps[[1]]
probes.near.genes <- probe.gene.overlaps[[2]]
probe.gene.map <- probe.gene.overlaps[[3]]

Fisher test against Gene Ontology and KEGG gene sets

Using topGO and topKEGG from limma.

Pathways enriched in genes close to any adipose DMRs

## Get enriched pathways etc. in genes close to hypermethylated adipose DMRs
enriched.terms.a.no.fc <- dmr.go.kegg(dmrcate.minfi.dmrs.nofc.a, close.to.probe.ensembl, gene.annot.tab, annoData, fisher.cutoff=0.01, min.genes = 3)

Pathways enriched in close to hypomethylated adipose DMRs

enriched.terms.a.down.no.fc <- dmr.go.kegg(dmrcate.minfi.dmrs.nofc.a[dmrcate.minfi.dmrs.nofc.a$meanbetafc <0,], 
                                                                close.to.probe.ensembl, gene.annot.tab, annoData, fisher.cutoff=0.01, min.genes = 3)

Pathways enriched in close to hypermethylated adipose DMRs

enriched.terms.a.up.no.fc <- dmr.go.kegg(dmrcate.minfi.dmrs.nofc.a[dmrcate.minfi.dmrs.nofc.a$meanbetafc >0,], 
                                                                close.to.probe.ensembl, gene.annot.tab, annoData, fisher.cutoff=0.01, min.genes = 3)

Enrichr

Enrichr ia a web tool where you can upload lists of genes, and look for overlaps with gene sets from many different sources (Gene Ontology, KEGG, Wikipathays, ENCODE, GTEx). The tool can be found at http://amp.pharm.mssm.edu/Enrichr/. The following code generates gene lists that can be pasted into the web form.

write.table(dmrcate.minfi.dmrs.a[dmrcate.minfi.dmrs.a$meanbetafc > 0,]$gene_tss, row.names = F, col.names = F, quote = F)

## TNXB
## ABAT
## PLEKHB1
## PHKG1
## LGI1
## PRDM1
## LGALS3BP
## NRXN1
## C8orf86
## RIN2
## PCDH15
## SOD3
## FZD4
## VMP1
## BAIAP2
## PRLR
## DCN
## RIMKLB
## IFLTD1
## CHD3
## FOXP2
## 
## TRIM2
## AKR1CL1
## ADAM5
## ASB16
## RP11-820L6.1
## 
## SPG20
## LINC00917
## RP11-209A2.1
## SULT1C2
## SEMA5B
## HMGB1P13
## CD36
## SLC38A4
## TEX26-AS1
## CTC-529P8.1
## 
## SMG5
## BCHE
## 
## SYNPO
## TMEM229B
## STARD8
## PLCD4
## CES3
## RNF183
## SAPCD1-AS1
## 
## SNAI2
## TMEM244
## COL11A2
## 
## CCDC80
## C7orf62
## HOXA2
## FXYD3
## PLA2G3
## CHM
## RGS12
## ST5

write.table(dmrcate.minfi.dmrs.a[dmrcate.minfi.dmrs.a$meanbetafc < 0,]$gene_tss, row.names = F, col.names = F, quote = F)

## TSPAN32
## 
## AVP
## WDFY4
## SNORD115-18
## ZNF385A
## TOP1MT
## RP11-627G18.2
## COL10A1
## ADORA2A
## UNCX
## MSX1
## 
## MAD1L1
## 
## CPT1A
## RHOD
## 
## RP11-441F2.5
## SLC25A30
## HOXD3
## RGS12
## TRIM39
## TMC6
## MPG
## RP11-109A6.3
## 
## HIC1
## GFI1
## MSX1
## CTSZ

write.table(dmrcate.minfi.dmrs.a$gene_tss, row.names = F, col.names = F, quote = F)

## TNXB
## TSPAN32
## 
## ABAT
## PLEKHB1
## PHKG1
## AVP
## LGI1
## WDFY4
## PRDM1
## LGALS3BP
## NRXN1
## C8orf86
## RIN2
## PCDH15
## SOD3
## FZD4
## VMP1
## SNORD115-18
## BAIAP2
## PRLR
## DCN
## ZNF385A
## RIMKLB
## IFLTD1
## CHD3
## FOXP2
## 
## TRIM2
## AKR1CL1
## ADAM5
## TOP1MT
## ASB16
## RP11-820L6.1
## 
## SPG20
## RP11-627G18.2
## LINC00917
## COL10A1
## RP11-209A2.1
## SULT1C2
## SEMA5B
## ADORA2A
## HMGB1P13
## CD36
## SLC38A4
## UNCX
## TEX26-AS1
## MSX1
## CTC-529P8.1
## 
## 
## MAD1L1
## 
## SMG5
## CPT1A
## RHOD
## 
## RP11-441F2.5
## SLC25A30
## BCHE
## 
## SYNPO
## HOXD3
## TMEM229B
## STARD8
## PLCD4
## CES3
## RNF183
## SAPCD1-AS1
## RGS12
## TRIM39
## 
## SNAI2
## TMEM244
## TMC6
## COL11A2
## 
## MPG
## CCDC80
## C7orf62
## HOXA2
## RP11-109A6.3
## FXYD3
## PLA2G3
## 
## HIC1
## GFI1
## MSX1
## CTSZ
## CHM
## RGS12
## ST5

Uploading the lists to Enrichr gave the following results:

• Adipose down http://amp.pharm.mssm.edu/Enrichr/enrich?dataset=3hvym

• Adipose up http://amp.pharm.mssm.edu/Enrichr/enrich?dataset=3hvxp

• Adipose all http://amp.pharm.mssm.edu/Enrichr/enrich?dataset=3hvz1

Print tables with DMRs and probes

write.table(as.data.frame(dmrcate.minfi.dmrs.a[dmrcate.minfi.dmrs.a$meanbetafc > 0,]), sep="\t", quote=F, row.names=F, file="./minfi_dmrcate_a_up.txt")
write.table(as.data.frame(dmrcate.minfi.dmrs.a[dmrcate.minfi.dmrs.a$meanbetafc < 0,]), sep="\t", quote=F, row.names=F, file="./minfi_dmrcate_a_down.txt")

write.table(as.data.frame(dmrcate.minfi.dmrs.nofc.a[dmrcate.minfi.dmrs.nofc.a$meanbetafc > 0,]), sep="\t", quote=F, row.names=F, file="./minfi_dmrcate_nofc_cutoff_a_up.txt")
write.table(as.data.frame(dmrcate.minfi.dmrs.nofc.a[dmrcate.minfi.dmrs.nofc.a$meanbetafc < 0,]), sep="\t", quote=F, row.names=F, file="./minfi_dmrcate_nofc_cutoff_a_down.txt")

Here we summarize the probe data by gene. For each gene, we use the mean M value of all probes annotated to that gene.

mvals.per.gene <- function(GRset, probe.gene.map, use.function=mean){
    filtered.map <- lapply(probe.gene.map, 
                 function(x){
                    intersect(x, rownames(GRset))
                 })
    
    mvals <- getM(GRset)
    out.list <- lapply(names(filtered.map),
                                        function(ensg){
                                            probe.data <- mvals[filtered.map[[ensg]], ]
                                            
                                            if(length(nrow(probe.data)) ==0){
                                                return(NULL)
                                            }
                                            
                                            return( apply(probe.data, 2, use.function, na.rm=TRUE) )
                                        })
    names(out.list) <- names(filtered.map) 
    out.tab <- do.call(rbind, out.list)
    return(out.tab)
    }
# This takes around 1 hour for each tissue, so it's commented out at the moment.
# mvals.per.gene.adipose <- mvals.per.gene(minfi.norm.a, probe.gene.map, use.function=mean)
# mvals.per.gene.muscle <- mvals.per.gene(minfi.norm.m, probe.gene.map, use.function=mean)
#
# write.csv(mvals.per.gene.adipose, file="gene_mean_mvals_adipose.csv",  quote=FALSE)
# write.csv(mvals.per.gene.muscle, file="gene_mean_mvals_muscle.csv",  quote=FALSE)

Finally, we export the M values as csv files. Since all probes overlppaoing SNPs have been filtered out, these data can be shared freely.

write.m.vals(minfi.norm.a, probes.near.genes, "./gene_mvals_adipose.csv.gz")
write.m.vals(minfi.norm.m, probes.near.genes, "./gene_mvals_muscle.csv.gz")

write.m.vals(minfi.norm.a, rownames(getM(minfi.norm.a)), "./non_snp_mvals_adipose.csv.gz")
write.m.vals(minfi.norm.m, rownames(getM(minfi.norm.m)), "./non_snp_mvals_muscle.csv.gz")

Session info

R.Version() 

## $platform
## [1] "x86_64-apple-darwin15.6.0"
## 
## $arch
## [1] "x86_64"
## 
## $os
## [1] "darwin15.6.0"
## 
## $system
## [1] "x86_64, darwin15.6.0"
## 
## $status
## [1] ""
## 
## $major
## [1] "3"
## 
## $minor
## [1] "4.3"
## 
## $year
## [1] "2017"
## 
## $month
## [1] "11"
## 
## $day
## [1] "30"
## 
## $`svn rev`
## [1] "73796"
## 
## $language
## [1] "R"
## 
## $version.string
## [1] "R version 3.4.3 (2017-11-30)"
## 
## $nickname
## [1] "Kite-Eating Tree"

sessionInfo(package = NULL)

## R version 3.4.3 (2017-11-30)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
##  [1] grid      splines   stats4    parallel  stats     graphics  grDevices
##  [8] utils     datasets  methods   base     
## 
## other attached packages:
##  [1] bindrcpp_0.2.2                                    
##  [2] gridExtra_2.3                                     
##  [3] pheatmap_1.0.8                                    
##  [4] biomaRt_2.34.2                                    
##  [5] EnsDb.Hsapiens.v75_2.99.0                         
##  [6] ensembldb_2.2.2                                   
##  [7] AnnotationFilter_1.2.0                            
##  [8] GenomicFeatures_1.30.3                            
##  [9] AnnotationDbi_1.40.0                              
## [10] ChIPpeakAnno_3.12.7                               
## [11] VennDiagram_1.6.20                                
## [12] futile.logger_1.4.3                               
## [13] RColorBrewer_1.1-2                                
## [14] DMRcate_1.14.0                                    
## [15] DMRcatedata_1.14.0                                
## [16] DSS_2.26.0                                        
## [17] bsseq_1.14.0                                      
## [18] limma_3.34.9                                      
## [19] doParallel_1.0.11                                 
## [20] sva_3.26.0                                        
## [21] BiocParallel_1.12.0                               
## [22] genefilter_1.60.0                                 
## [23] mgcv_1.8-23                                       
## [24] nlme_3.1-137                                      
## [25] minfiData_0.24.0                                  
## [26] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.0
## [27] IlluminaHumanMethylation450kmanifest_0.4.0        
## [28] minfi_1.24.0                                      
## [29] bumphunter_1.20.0                                 
## [30] locfit_1.5-9.1                                    
## [31] iterators_1.0.9                                   
## [32] foreach_1.4.4                                     
## [33] Biostrings_2.46.0                                 
## [34] XVector_0.18.0                                    
## [35] SummarizedExperiment_1.8.1                        
## [36] DelayedArray_0.4.1                                
## [37] matrixStats_0.53.1                                
## [38] Biobase_2.38.0                                    
## [39] GenomicRanges_1.30.3                              
## [40] GenomeInfoDb_1.14.0                               
## [41] IRanges_2.12.0                                    
## [42] S4Vectors_0.16.0                                  
## [43] BiocGenerics_0.24.0                               
## [44] forcats_0.3.0                                     
## [45] stringr_1.3.1                                     
## [46] dplyr_0.7.5                                       
## [47] purrr_0.2.4                                       
## [48] readr_1.1.1                                       
## [49] tidyr_0.8.1                                       
## [50] tibble_1.4.2                                      
## [51] ggplot2_2.2.1                                     
## [52] tidyverse_1.2.1                                   
## [53] data.table_1.10.4-3                               
## [54] knitr_1.20                                        
## 
## loaded via a namespace (and not attached):
##   [1] R.utils_2.6.0                                      
##   [2] tidyselect_0.2.4                                   
##   [3] IlluminaHumanMethylationEPICanno.ilm10b2.hg19_0.6.0
##   [4] RSQLite_2.1.1                                      
##   [5] htmlwidgets_1.2                                    
##   [6] munsell_0.4.3                                      
##   [7] codetools_0.2-15                                   
##   [8] preprocessCore_1.40.0                              
##   [9] statmod_1.4.30                                     
##  [10] colorspace_1.3-2                                   
##  [11] BiocInstaller_1.28.0                               
##  [12] rstudioapi_0.7                                     
##  [13] labeling_0.3                                       
##  [14] GenomeInfoDbData_1.0.0                             
##  [15] mnormt_1.5-5                                       
##  [16] bit64_0.9-7                                        
##  [17] rprojroot_1.3-2                                    
##  [18] lambda.r_1.2.2                                     
##  [19] biovizBase_1.26.0                                  
##  [20] regioneR_1.10.0                                    
##  [21] R6_2.2.2                                           
##  [22] illuminaio_0.20.0                                  
##  [23] idr_1.2                                            
##  [24] bitops_1.0-6                                       
##  [25] reshape_0.8.7                                      
##  [26] assertthat_0.2.0                                   
##  [27] promises_1.0.1                                     
##  [28] scales_0.5.0                                       
##  [29] nnet_7.3-12                                        
##  [30] gtable_0.2.0                                       
##  [31] methylumi_2.24.1                                   
##  [32] rlang_0.2.0                                        
##  [33] rtracklayer_1.38.3                                 
##  [34] lazyeval_0.2.1                                     
##  [35] acepack_1.4.1                                      
##  [36] GEOquery_2.46.15                                   
##  [37] dichromat_2.0-0                                    
##  [38] broom_0.4.4                                        
##  [39] checkmate_1.8.5                                    
##  [40] yaml_2.1.19                                        
##  [41] reshape2_1.4.3                                     
##  [42] modelr_0.1.2                                       
##  [43] backports_1.1.2                                    
##  [44] httpuv_1.4.3                                       
##  [45] Hmisc_4.1-1                                        
##  [46] RBGL_1.54.0                                        
##  [47] RMySQL_0.10.15                                     
##  [48] tools_3.4.3                                        
##  [49] psych_1.8.4                                        
##  [50] nor1mix_1.2-3                                      
##  [51] siggenes_1.52.0                                    
##  [52] Rcpp_0.12.16                                       
##  [53] plyr_1.8.4                                         
##  [54] base64enc_0.1-3                                    
##  [55] progress_1.1.2                                     
##  [56] zlibbioc_1.24.0                                    
##  [57] BiasedUrn_1.07                                     
##  [58] RCurl_1.95-4.10                                    
##  [59] prettyunits_1.0.2                                  
##  [60] rpart_4.1-13                                       
##  [61] openssl_1.0.1                                      
##  [62] haven_1.1.1                                        
##  [63] cluster_2.0.7-1                                    
##  [64] magrittr_1.5                                       
##  [65] futile.options_1.0.1                               
##  [66] ProtGenerics_1.10.0                                
##  [67] missMethyl_1.12.0                                  
##  [68] hms_0.4.2                                          
##  [69] mime_0.5                                           
##  [70] evaluate_0.10.1                                    
##  [71] xtable_1.8-2                                       
##  [72] XML_3.98-1.11                                      
##  [73] mclust_5.4                                         
##  [74] readxl_1.1.0                                       
##  [75] compiler_3.4.3                                     
##  [76] crayon_1.3.4                                       
##  [77] R.oo_1.22.0                                        
##  [78] htmltools_0.3.6                                    
##  [79] later_0.7.2                                        
##  [80] Formula_1.2-3                                      
##  [81] lubridate_1.7.4                                    
##  [82] DBI_1.0.0                                          
##  [83] formatR_1.5                                        
##  [84] MASS_7.3-50                                        
##  [85] ade4_1.7-11                                        
##  [86] Matrix_1.2-14                                      
##  [87] permute_0.9-4                                      
##  [88] cli_1.0.0                                          
##  [89] quadprog_1.5-5                                     
##  [90] R.methodsS3_1.7.1                                  
##  [91] Gviz_1.22.3                                        
##  [92] bindr_0.1.1                                        
##  [93] pkgconfig_2.0.1                                    
##  [94] GenomicAlignments_1.14.2                           
##  [95] registry_0.5                                       
##  [96] foreign_0.8-70                                     
##  [97] xml2_1.2.0                                         
##  [98] annotate_1.56.2                                    
##  [99] rngtools_1.3.1                                     
## [100] ruv_0.9.7                                          
## [101] pkgmaker_0.22                                      
## [102] multtest_2.34.0                                    
## [103] beanplot_1.2                                       
## [104] rvest_0.3.2                                        
## [105] doRNG_1.6.6                                        
## [106] VariantAnnotation_1.24.5                           
## [107] digest_0.6.15                                      
## [108] graph_1.56.0                                       
## [109] rmarkdown_1.9                                      
## [110] base64_2.0                                         
## [111] cellranger_1.1.0                                   
## [112] htmlTable_1.11.2                                   
## [113] curl_3.2                                           
## [114] shiny_1.0.5                                        
## [115] Rsamtools_1.30.0                                   
## [116] gtools_3.5.0                                       
## [117] jsonlite_1.5                                       
## [118] seqinr_3.4-5                                       
## [119] BSgenome_1.46.0                                    
## [120] pillar_1.2.2                                       
## [121] lattice_0.20-35                                    
## [122] GO.db_3.5.0                                        
## [123] httr_1.3.1                                         
## [124] survival_2.42-3                                    
## [125] interactiveDisplayBase_1.16.0                      
## [126] glue_1.2.0                                         
## [127] bit_1.1-13                                         
## [128] stringi_1.2.2                                      
## [129] blob_1.1.1                                         
## [130] org.Hs.eg.db_3.5.0                                 
## [131] AnnotationHub_2.10.1                               
## [132] latticeExtra_0.6-28                                
## [133] memoise_1.1.0